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Introduction 



In doing numerical weather forecasting there are essentially two 
problems: 1) the forecasting, which amounts to solving a coupled set of 

non-linear partial differential equations in space and time, and 2) 
determining the initial conditions for those equations from observational 
data. The problem of determining the initial conditions accurately is a 
very serious one which has not yet been approached by finite element 
methods. Thus we shall not discuss it in this report. 

The forecast equations are essentially the laws of conservations of 
mass, momentimi and energy. Written in the usual form, as a set of first 
order partial differential equations, they are called the "Primitive 
Equations" or the "P.E. model." As such they are hyperbolic, and non- 
linear. These equations are 



(tT + V*V + w |-)V + FkxV + -Vp= vV^V (la) 

dt dZ p 

+ + w |^)log 0 = (Ic) 

•^ + V* (pV) + (pw) = 0 (Id) 

3t dz 



where F is the Coriolis "parameter" 2 sin 0, 0 is the potential 

-R/Cp 

temperature (0 = T(P/P^) ), V is the horizontal velocity vector 

and w is the vertical velocity. The wave solutions which result 
consist of two types of waves - the so-called inertial-gravity waves and 
the Rossby waves (See [1]). It is known that the gravity waves should 
not be present, and in any numerical computation they die out, with time 



- 1 - 



constants of 12-24 hours. However, this length of time is the time of 



primary interest for forecast purposes. 

The laws of conservation can also be written in terms of scalar and 
vector potentials i.e., a geopotential 4^, a vorticity 5 and a 
stream function }p . These functions satisfy a second order elliptic 
and a single hyperbolic equation and thus have the property that the 
gravity waves are damped out, or filtered. (See equations (27) and (29)). 
This form of the equations is referred to as the vorticity, or filtered, 
form of the equations. However, involving a second order equation they 
require different boundary conditions, which conditions are somewhat 
contrived. 

For the primitive equations (1), one can set up an order of magnitude 
analysis and determine which terms are significant for meteorological 
purposes. This leads to a hierarchy of equations which one can solve 
more easily than (1). The zero order analysis gives the following set 



which dominate all the other terms by an order of magnitude. This 
suggests that the basic flow is horizontal, divergence free, and driven 
by the balance of the pressure force and the Coriolis force. 

At the next level we have the so called barotopic equations 



F k X V = - 




(2a) 



V-(pV) = 0, 



(2b) 



(|- + V-V)V + FkxV+V<}. = 0 

a t 



(3a) 



+ 4 >v*v = 0 



(3b) 
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where (j) = gh is the geopotential at the height h, obtained by inte- 
grating dp in the z direction, since equations (3) are independent 

of vertical variables. 

If we take a perturbation approach about some uniform velocity 
and a height H, we get a set of hyperbolic first order equations 



|^=-V *Vu + Fv-|^ 
3t o 9x 



‘Vv- Fu-|^ 
9t o 9y 



li - gH (iii + 

9t o ^ ^ o^9x 9y-^ 



A A 

where V = iu + jv . If we look for a solution of the form 



\p = (u,v,(f>) = Ijj e 
o 



i (wt-k*x) 



we find three waves 



oj = V • k 
o o 



^ f Ak A. 

• k ±VgH k*k + F 



“1,2 = ^0* ^ ±^®«o 



(4a) 

(4b) 

(4c) 



(5) 



(6a) 

(6b) 



The first is a convective wave, due to the uniform flow. The second and 



^ li' 

3x 



) were zero 



third are **inertial-gravity** waves. If the terms (Fv, 

in (4a) and (-Fu, were zero in (4b), then these twc waves would 

oy 

disappear. These are just the terms in (2a), which essentially "balance 



out. 



If the initial conditions contain any of the eigenvectors to which 
these last two waves are the eigenvalues, then the linear solution 
contains them. Thus any solution to the non-linear problem which 
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initially contains some component in these ^directions** will propagate 
them for some time. The real problem is non-linear, and in actuality 
they die out but we would hope that they are not there to start with. 

As a result there is considerable interest in getting correct 
initial conditions from the observational data. We turn first to the 
methods which involve vorticity. These all involve the assumption that 
there is a static balance of the wind and pressure fields, and that one 
can be derived from the other. Proceeding from there is complicated 
by the fact that in the tropics there are very small pressure changes, 
so that observational errors in pressure can be as great as the pressure 
changes themselves, while the winds tend to be more accurately recorded. 
The converse is true in the mid - and high - latitudes. So in the 
raid-latitudes one takes the pressure measurements at various points 
(invariably not grid points), computes the pressure at the grid points 
by some weighted average, converts to the geopotential cf), and then 
computes the stream function using one of the three relations (9) - 
(11). The winds are then obtained from In the tropics the 

proceedure is reversed. The observed winds V are used to compute the 
vorticity 



Then the geopotential ({> is obtained by solving one of the following: 



C = k-V X V 



(7) 



The stream function ip is obtained from the Poisson equation, 




( 8 ) 



2 2 
V^(f) = FV ip 



(9) 
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( 10 ) 




( 11 ) 



where F is again the Coriolis factor. 



T - 2Q sin0 



( 12 ) 



J is the symbol for the Jacobian, and 0 is the latitude from the 
equator. 

These equations are called the quasi-geostrophic relation, the 
linear balance equation, and the non-linear balance equation. Note from 
(12) that to use (9), (10) or (11) to solve for at the equator 
involves solving a singular differential operator. Note also that these 
equations are all Poisson equations. Note also that this procedure only 
deals with one part of the wind, the rotational part. Nowhere does the 
divergent part enter the calculation. Since the wind is assumed 
non-divergent in the first place, that is not essential at the start. 
However, a finite difference scheme will introduce some errors, and this 
may be significant. 

2. Finite Difference Models 

These methods were used, and are still used in some places, for 
many years. However, more and more people are turning to the primitive 
equations. So let us consider the simplest version, the barotropic 
model. First we note that the gravity waves (6b) have a much higher 
speed (about 300 meters/second) than the meteorological mode (6a) (5 
meters/second). For computing , using a system of hyperbolic equations, 
the Courant-Lewy-Fredericks criterion requires that the time step be 
less than Ax divided by the wave speed. These fast waves then dominate 
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the size of the time step, for calculating using an explicit calculation. 
Thus there is considerable interest in semi-implicit schemes. 

Williamson and Browning [2] run some computer tests on the primitive 
equations, using the barotropic model* writing the equations in spherical 
coordinates. They use both the advection form, (a is the Earth's radius^ 



9u u 9u V 9u , . u . 1 9d) 

T T:; — ~ TT + (f i tan 0)v 

9t acos0 9X a 90 a acos0 9X 



9v u 9v V 9v /^ . u ^ 1 94> 

T" = - T (f+“ tan 0)u — ^ 

9t acos0 9X a 90 a a 90 



{•^(<J)U) + -^(4>V COS0)} 



9t 



acos0 9X 



90 



(13) 

(14) 

(15) 



and the flux form of (13) and (14) 



+ a tan0)<^v 






(16) 



+ |^(<^v^cos0)} - (f + 7 tan0)<^u 

(17) 



In their numerical scheme they use a centered time difference and a 

centered space difference for all derivatives. For the terms involving 

(f + — tan0) they use a time average (u = -^ [u(t+l, x) + u(t-l, x)]) 
a z 

so that these terms in either (13) and (14) or (16) and (17) are treated 
implicitly - thus creating a semi-implicit model. 

They ran three tests. In the first, they used 5° grids in both X 
and 0, taking the time step based on the smallest spacial distance near 
the pole. In the second they arbitrarily omitted points on longitude 
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circles as the latitude increased towards the pole, so that at 60°N they 
had 72 points per longitude circle and at 85°N they had 12. (This has 
the same effect as an icosohedral grid). This allows a longer time step, 
by a factor of 6. In the third, they took the output, did a longitudinal 
Fourier analysis and threw away the high frequency, short wave components 
("filtering"). (This is what is proposed for FNWC P.E. model). They 
start up the scheme with a knoxm steady state solution of (13) - (15), 
namely zonal geostrophic flow, given by 

u = u (cos0 cosa - cosX sin0 sina) 
o 

V = u sinX sina ^ (18) 

o 2 ^ ^ 

^o 2 

(t> = (p^ - (a Q -1 ^) (cosX COS0 sina + sin0 cosa) 

(If a = 7 t/2 there is no tendency to flow across the equator, which 
simplifies the calculation). 

Their results are - for test one 

1) For the initial time step the error in u is dominated by the 
truncation error in (f>. 

2) The error in u grew by two orders of magnitude in five days, 
from the error after one time step. 

3) The errors start at the pole and spread out. 

In test two, there was 

1) a factor of ten in the size of the error near the poles. 

In test three they used the same time step as in test two (six minutes 
vs one minute for test one). The computations were three times faster 
than in test one. They were able to reproduce exactly the results of 



test one. With a fourth order difference scheme in space and the filtered 



scheme they were able to get the best results by about two orders of 
magnitude. Thus for this type calculation a filtered model is clearly 
the best. 

Kwizak and Roberts [3] rewrite the three equations (13) - (15) as 
follows: 

Let 



and 



Then we can write 



1 . 2 
2 



^ * lie 'i 



Ju 

8t 



1 

acos0 9X 



at 



a ae 



at 



1 r9 



acos6 ax 



+ v^) 


(19) 


- — (u cose)]. 


(20) 


0 1 3K 

acos0 ax 


(21) 




(22) 


+ ■|^('))V COS0)} 


(23) 



The function K is interpretable as the kinetic energy and Q is 
interpretable as the absolute vorticity [3] or the potential vorticity 
[4]. The latter comes from the fact that (20) can be written 



Q = F + V X V (24) 

where V is the velocity vector. 

Now the original equations are derivable from assuming an 
incompressible gas, so that we have the implicit equation 
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which implies that there is a stream function i}) such that 



u 



dll; 

dy ^ 



V 



ii 

9x 



so that (23) can be rewritten 



V 




Q - 



F 



And taking the curl of (21) - (22) we get 



- u dQ V dQ ^ 

9t ^ acos9 9X ^ a 30 ° 



which can be rewritten 



= J(Q,\l>) 



( 26 ) 



(27) 



(28) 



(29) 



where J is the Jacobian of Q and ip. The pair (26) and (29) are a 

coupled system, with (26) being an elliptic equation, which could be 

solved iteratively, and then (j) could be found from (15). 

Observe that taking the curl of (21) - (22) eliminates the <t> from 

these equations. Typical values of the variables in the atmosphere are 
4 

(|)-3xl0,u=5,v=l. So a potentially large term has been 

eliminated from the computations for the velocity V. Also observe that 

(28) is the only wave-like equation left for V and this, being first 

order, has only one wave solution. A linearization of (28) about a 

steady flow will give only one wave 

ik(x-V t) 

^ o 



which is the ^’meteorological” wave which is expected: that is, the two 

gravity waves do not enter into the computation for Q. This is the 
reason that the vorticity form of the equations was used for many years - 
the time dependent portion of the problem involved only a single wave 
speed, the one which was desired. The problem is that the solution must 
be kept divergence free; this may be a nontrivial task. For experimental 
purposes we can see that the ’’meteorological wave” ought to be associated 
with having a divergence free field, and thus making an effort to get the 
initial conditions divergence-free is a worthwhile computation, as far 
as reducing ’’unwanted noise” in the answer. For this reason Kwizak and 
Roberts comment ”... the winds are perfectly non-divergent' initially and 
at the end of the first time step. This property virtually eliminates 
the gravity waves from the integration,” [3]. 

The computations which Kwizak and Roberts do are based on the three 
primitive equations (21), (22) and (15), doing a semi-implicit scheme. 
They take the (|) terms in (21) and (22), and then rewrite (15) as 

and then treat the first term implicitly. They thus have three coupled 

equations to handle. The actual method of computation is to observe 

that (21) and (22) can be used to convert (29) to a non-homogeneous 

-2t -2t 

Helmholtz equation for (j) , defined as ^ = [(|)(t+At) + (j)(t-At)], 

namely 

= f(*.u.v) (31) 
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Having found this, then (21) and (22) can be used to compute u and v 
in a purely explicit fashion. For their computations they can use a 60 
minute time step (vs a 10 minute for an explicit calculation) at a 
savings of 4 times in the computation speed. Thus the method of 
computation does not involve three coupled implicit equations but one 
Helmholtz equation (second order) and two explicit calculations. 

Two years later Elviers and Sundstrom [5] do a similar test, using 
essentially the same scheme, noting that the averaging operator which is 
used allows for a decoupling of the equations for even and odd time steps, 
i.e., a staggered grid. They do a stability analysis of the finite 
difference models. Their analysis shows the semi-implicit scheme is 
superior to an explicit scheme. 

Williamson and Browning [2] and Kwizak and Roberts [3] both do a 
semi- implicit scheme, but they differ as to which terms they treat 
implicitly. Kwizak and Roberts do the more usual method. In the u and 
V equations they treat the pressure gradient term (involving <()) 
implicitly. This allows them to eventually eliminate u and v from 
the (j) equation and convert the latter to a second order Helmholtz 
equation. This one equation is solved fully implicitly, with u and v 
then found explicitly. Williamson and Browning treat the Coriolis terms 
implicitly, giving a coupled system for u and v, while the ^ 
equation is treated explicitly. In both cases only one of the terms 
driving the gravity waves is treated implicitly, but this is sufficient 
to remove them from the Courant criterion. 

Thus a significant problem for the predictive equations is the 
ability to solve quickly and efficiently a system of partial differential 
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equations. The usual method is to use classical finite difference 
techniques. The standard way is to use second order centered differences. 
There is every reason to believe that fourth order centered differences in 
space will increase the accuracy. This is currently being investigated [6]. 

3. Finite Element Methods 

A question of definite interest is whether the "finite element 
method" of solving partial differential equations will give better 
accuracy. I am aware of only four places where this is currently being 
investigated. George Fix [4,7] is studying ocean circulation problems 
this way. His studies are being continued by Hirsch [17]. M. P. Cullen 
[8, 9, 13] has programmed the barotropic equations, and is now attempting 
to program a more realistic set of primitive equations. A. Staniforth 
[10] is attempting to implement finite element calculations in the 
Canadian Meteorological Office. And a student thesis at NPS by Donald 
Kinsman [11] ran some experiments with finite elements on the barotropic 
equations. All of these indicate that this method may have a significant 
future in meteorology. 

Fix [4] looks at the ocean circulation problem, which is just (13) 
and (14) , together with the divergence free condition 

V • V = 0 

Thus he does not get involved with gravity waves and has only the 
"meteorological mode" to contend with. He then converts to the vorticity 
equation and the non-linear balance equation (27), (29). He then takes a 
finite element - Galerkin approach, using linear elements (also 
quadratics and cubics for further tests). 
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There are three problems to be addressed in any analysis of this 

discretization. The first is the accuracy of the special discretization. 

The second, which is due to the fact that these equations are non-linear, 

is called "aliasing," a feature which can be most easily seen by 

considering a Fourier analysis of the special terms. If both u and v 

til 

are written as Fourier series and then truncated at the N— ^ harmonic, 

then a term involving u times v will have a term, say 

(u^cos N x) (v^cos M x), which would normally give rise to two terms 

^ ^ latter can not appear 

due to the truncation. Certain discretization schemes have an imbalance 

in the treatment of this phenomenon, knoxm as aliasing. However, 

Jesperson [12] has shown that this phenomenon does not occur with a 

finite element scheme, that is, finite element schemes are free from 

aliasing, a fact which Fix reconfirms. 

Fix also shows what is widely known, that the special accuracy for 
k 

the velocity is 0(h) where h is the grid size and k is 1, 2, or 
three depending upon whether linear, quadratic or cubic elements are 
used. 

The third problem is how to handle the time integration. Fix does 
not have to worry about a semi-implicit scheme from the point of view of 
gravity waves. However, any finite element scheme links more than two 
points and one is automatically forced into an implicit scheme. That is 
one drawback for finite elements. Fix proceeds to analyze a linear one 
dimensional analogue of the wave propogation problem. The finite 
element method generates its own natural set of difference equations. 

For the linear problem which Fix sets up, using linear elements in space 
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with time varying coefficients, Fix shows that the implicit discretization 
which the finite element method forces gives fourth order accuracy in the 
phase speed of waves, a fact which Cullen had noted earlier [8]« 

Fix then chooses for his time discretization the usual centered 
difference (leap-frog) time discretization, so that the final computation 
for Q is given by 



Thus for the non-linear ocean circulation problem one knows that 
the finite element method does not introduce aliasing, is spatially as 
accurate as we make the finite elements, and can be conjectured to be 
fourth order accurate in phase speed, (Phase speed of the waves has 
always been a problem in forecasting of weather). 

In a series of three papers [8], [9], [13], Cullen tackles two 
problems. The first problem is to solve equations (3) in a limited area 
-1 _< X, y ± 1 (the **beta plane”) with periodic boundary conditions, 
using a finite element method. The second is to solve equations (3) on 
the surface of the globe. His analysis proceeds as follows. 

He first [13] considers a single linear equation 



on a rectangle, using rectangles and bilinear elements, with V known. 

He compares his results, using a 16 x 16 grid, to second and fourth 
order finite difference schemes using a 32 x 32 grid. The exact solution 
to this problem is known. He demonstrates that the finite element 




A 



A 



+ V-V(^ = 0 



(32) 
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calculation is more accurate, both in handling smooth data and in handling 
a problem with a discontinuity. 

Satisfied with the results of the linear problem, Cullen [13] 
proceeds to the non-linear problem (3), in a grid -1 _< x, y Ji 1> without 
the Coriolis term. His initial condition would give a gravity wave if 
the problem were linear. He compares his answer with a 16 x 16 grid to 
finite difference schemes with 32 x 32 and 64 x 64 grids. He uses linear 
finite elements with a leap-frog time step. His results indicated that 
the finite element method was better, although there is no exact solution 
to compare with. Cullen then attempts to analyze the numerical scheme. 

For the linear one dimensional case he shows that the phase error is 
fourth order in At/Ax. In fact, he shows something more, namely that 
for some range of At/Ax the phase error can actually cause a small 
leading phase. He also argues, largely on qualitative grounds, that 
there is no aliasing problem. There is no discussion of the inertial vs 
Rossby waves, and no attempt to isolate one from the other, or to control 
either. 

Cullen’s second paper [8] concerns equations (3) on the entire 

globe - the genuine barotropic problem. One question which Cullen now 

addresses is how is the best way to handle the non-linear terms which 

appear in (3). His approach is to analyze a single one dimensional term 

of the form u . The various possibilities are: 
ox 

1. Treat v as known, take the nodel values of v, compute a 
derivative, and simply multiply the coefficients for u by these values. 

2. Compute a finite element expansion for v, analytically 

3v 

differentiate this expansion to find — , and use the resultant two 

ox 

expansions to compute t = u . 
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3. Compute a finite element expansion for — itself. 

Cullen claims that the last method, while sacrificing some accuracy, 
controls aliasing much better than either of the first two. 

I was unable to follow his arguments, and will attempt in a further 
report to see whether it can be generalized. In any event Cullen uses 
the third method in his calculations [8]. 

In [8] Cullen runs four tests on the equations (3), using the three 
methods above (for two of the runs he slightly modifies the coefficients 
in method 3) . He takes as initial values a finite element solution to 
(2a), comparing with other published results using finite differences. 
Assuming that the published finite difference results are the most 
accurate (?) he notes that the phases on most finite element runs appear 
to lag, although some are advanced. He concludes that: 

(a) Waves down to four element lengths will be treated almost 
perfectly. 

(b) Waves less than two grid lengths will not be treated at all well. 

(c) The finite element scheme ’’essentially eliminates aliasing 
errors . ” 

(d) The boundary condtions used on the problem introduced errors of 
the same order of magnitude as the change from finite differences to 
finite elements, and are thus quite significant. 

In his third paper [9] Cullen reports on actual computations on a 
sphere, relying heavily on the analysis above. He concludes ’’the finite 
element method is computationally more efficient.” He uses an icosahedral 
grid, subdivided by latitude and longitude lines to form triangles, 
resulting in 1002 points. In integration he treats the trig functions in 
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(3) as constants over each small triangle. He uses the scheme of [8] to 
treat the non-linear terms. 

For his actual computational scheme it appears that the only implicit 
portion of the scheme is the implicitness generated by the left hand side. 
In other words, looking at (3) as equations of the form 

f - ’'(♦> 

he solves the resultant non-linear equations by the following iterative 
process. Treat the right hand side as known and take the diagonally 
dominant matrix on the left as the generator of the next iteration. In 
this way the iterative process to find the values of the nodal points at 
time t is relatively fast. Then a leap-frog time discretization of 10 
minutes was used. (He could use one of up to 14 minutes). He found that 
filtering was required every two hours to get long time (greater than 5 
days) solutions. His initial conditions were Rossby waves with wave 
number 4 and wave number 8. He compares his results with published 
results using a finite difference model with 4032 points, one with 14,592 
points and a spectral model with 640 degrees of freedom. His results are 
better than the first but not as good as the last two. He also observes 
that the errors seem to start from the vertices of the icosahedron, where 
there are only 5 supporting triangles instead of 6 as there are at the 
intermediate triangulation points. 

Hinsman [11] in his master's thesis again studied equations (3). He 
considered two possible grids. One used the two angles A and 6 as 
rectangular coordinates and triangulating the resultant "rectangle." 

This results in a very fine subdivision of the polar regions resulting in 
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very small physical lengths Ax in that region. The second grid was an 
icosahedral grid, (as did Cullen [8], [9]) which was subsequently sub- 
divided by arcs of great circles. This results in most of the nodes 
having 6 adjacent nodes. The corners of the icosahedron have only 5 
adjacent nodes. This appears to generate some **noise,” as Cullen noted 
[9]. 

Experiments with the (X,0) grid showed instability after 12 hours, 
as might be guessed from the fact that there were 36 points around each 
latitude circle including those of the poles. The instability clearly 
arose in the polar region. Experiments with the icosahedral grid did not 
show these problems. 

Starting with an analytic solution to the non-linear balance equation 
which has essentially one wave going around the earth. Kinsman Fourier 
analyzed the solution as it propagated, comparing his results to similar 
finite difference calculations. In the low latitudes the propagation 
speed was almost exactly correct; as compared to 50-60% correct for 
finite differences. In the high latitudes both finite elements and 
finite differences fall behind the predicted speed. (This contradicts 
Cullen’s observations) . 

The method of solution was very different from any previous methods. 
Each equation in (3) is ”quasi-linearized” by considering the other two 
variables as known (from a previous calculation) . The finite element 
scheme automatically generates an implicit scheme, but the decision was 
made to go one step further and use a Crank-Nicholson approach, 
essentially calculating the variables at time step (N+1/2). The quasi- 
linearization uncoupled the equations, but, like Cullen, there was no 
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attempt to distinguish between the Rossby wave and the gravity waves. 
Finite elements were used to expand all the variables, having first 
integrated by parts to get the weak form. The equations were then 
successively solved, first for <(>, then for u and finally for v. 

(The trig functions sin0 and cos0 were also written as finite 
elements, as opposed to Cullen who treated them as constants over 
triangles). The resultant equations were solved by a Gauss-Seidel 
iteration. This took about 10-12 iterations to converge. (An additional 
feature of the program was a very efficient coding scheme to avoid 
completely the storage of the zero elements of the matrices involved). 

Swartz and Wendroff [14] compare the relative efficiency of finite 

difference and finite element methods. They do so by taking a one 

dimensional linear problem. Their interest is somewhat different than 

ours and thus their data is recorded very differently. They record, for 

example, (Table 1) the number of intervals per wave length which are 

necessary to get a desired accuracy in phase (in theory). One conclusion 

which can be drawn is that a dramatic improvement in resolution can be 

obtained by switching from linear elements to quadratic elements (for an 
-2 

error of 10 the number of intervals per wavelength goes from 8.7 to 
4.8, while for 10 ^ error the change is from 27 to 9.7). 

They also, using results of Kreiss and Oliger, report the results of 
finite difference theory. For fourth order special accuracy, - which is 
what linear elements give - the corresponding results are a change from 
13.3 to 7.9 and from 42.5 to 17.3. So one can conclude that in theory 
the finite elements have better phase resolution (e.g., 8.7 intervals per 
wave length vs 13.3 for linear elements and 4th order differencing) and 
that the payoff for increasee complexity is also greater. 
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The method that they use for actually computing the finite element 
solutions is interesting - they call it ”use of the trapezoidal rule.” 
Consider the equation 

d\x _ du 
3t ^ 8 x 



and compute the finite special element expansion of this to get a formula 

A u = B u (33) 

where A and B are sparse matrices and u is a vector (for linear 

elements A has a -7 (1, 4, 1) structure). Now take a time discretization 

6 

which gives 



A 



"'n+1 ^n 
u - u 



At 



'“n+l , ^n 
u + u 



(34) 



This scheme is fully implicit and is stable for all h, At. They compare 

this with the leap-frog finite difference scheme. 

Before finishing with the linear problem they conclude that a 

twelfth order special difference scheme, with leapfrog, is competitive 

— 4 

with a cubic spline, for phase errors of 10 

They turn to the first order one-dimensional non-linear problem 

Ut = f(u) + g(u) . (35) 

They have two conclusions, one of which is given just in passing. They 

remark that they will not use finite elements directly on (41). Instead 

—1 9 

they will use the operator A B to replace — , where A and B are 

dX 

defined in (33). Their reasoning is that the finite element scheme is 
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awkward and that it has been shown that use of the finite element scheme 



directly will degrade the special accuracy [15]* The scheme above will 
keep the truncation error. They further remark that **It no longer makes 
sense to try to find the best scheme. A reasonable approach seems to be 
to compare schemes using the same number of intervals per wave length." 
This sort of reasoning is directly counter to what Cullen and Kinsman 
have used. 

Their second conclusion, based on the above and counting evaluations 
of f(u) and g(u) as the major source of computer time., is that, using 
cubic splines and 18th order differencing as comparable schemes, finite 
difference and finite element schemes are roughly comparable. (For 
global meteorology an 18th order difference scheme seems excessive). 

Another paper which deals with non-linear finite element calculations 
for hyperbolic problems is by Oden and Post. What Oden and Post do to 
find the finite element formula for a hyperbolic problem is to use as a 



the analogue of — = 0 , where E is the energy of the system. It is 
at 

then clear that their method preserves energy. In the actual Computation 



test function V 



i 




so that the equation with which they work is 



dE 



9Vi 

they use the centered (leapfrog) approximation to -r— , and their 

o t 



analysis is of a single second order equation 




(36) 



which, in weak form, is 




0 
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Their study reveals the following: 

tl 

(a) the resultant equation is numerically stable if — ^(max)’ 



where C, s is the maximum wave speed C/T' (u ) . 
(max) X 

(b) If instead of using the matrix generated by 



2 ’ Iri 



so-called consistent mass matrix) one uses a diagonal matrix (the lumped 
mass approximation) the stability criterion is improved to read 

At ^(raax)- 

(c) The finite element solution converges uniformly to the solution, 
using piecewise linear elements. They do not make any analysis of the 
order of accuracy of this approach, so this paper does not touch the 
question raised by Swartz and Wendroff. 



4. Questions and Some Answers 

A number of questions arise from the intersection of these papers. 

1. What is the best method of handling the non-linear terms, and 
the variable coefficients? Cullen, Kinsman, Oden and Swartz all advocate 
different answers for different reasons. 

2. What does the finite element method do to the gravity waves 
which worry the people doing finite difference calculations? None of the 
finite element calculations even mention this question. 

3. What are the merits of the **lumped mass” vs ’’consistent mass” 
approach which is so familiar to the mechanical engineers? 

4. If this were to be implemented as an operational scheme, what 
are the possible merits of SOR or ADI. Hirsch [17], doing ocean 
circulation problems, solves the Poisson equation (27) by SOR and the 
advection equation (29) by ADI. Staniforth and Mitchell [10] use an ADI 
method for their calculations. 
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5. Cullen does not integrate by parts to compute his u — 

oX 

term. What effect does this have on the solution? Kinsman (unpublished) 
noted a considerable improvement in his results when he used the weak 
form (integrated by parts) of the (|) equation as opposed to the strong 
form (3). 

6. Fix has proved that the vorticity equations for ocean circulation 
automatically satisfy the desired conservation laws when written as 
finite elements. Is this also true for the barotropic equations? 

7. What does the finite element method do to the phase speed in 
two dimensions? 

We now proceed to answer this last question; at least for a linear 
model. The finite element approximation to 

+ V*V(j) = 0 (37) 

using bilinear elements on rectangles is 



kh 



35<16 



3hu 



‘*’)l+l,j-l ‘*’)l-l,j-l^^ 36 ‘*’^-+1,3 ” ^ ^‘*’«.+l,j+l ~ ‘**)l-l,j+l^ 



3kv 









We investigate the phase propagation of this linear model. The exponen- 
tial solution to (37) is ()) = A e ^ where oj = V*k. If we look 

* 

for (0 in the approximation 
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(38) 



<!> 






* 



ic 

^ ^-ico t + i£k+ijh 



then 



-ikho) {16 + 8 [cos h + cos k] + 4 cos k cos h} 

+ 3hui(sin k}(2 + cos h) + 3kvi{sin h}(2 + cos k) = 0 



thus 



* _ 3uh sin k(2 + cos h) + 3vk sin h(2 + cos k) 
^ hk(2 + cos h) (2 + cos k) 



3 sin k 
(2 + cos k)k 



3 sin h 

^ (2 + cos h)h 



(39) 



We then have the following theorem. 

Theorem 1: The finite element approximation to (37) using bilinear 

elements on rectangles has fourth order accuracy in phase speed. 

Proof: From (39) 



* 

0 ) - u 



( 3 - 1 + 31 , - ... ) 

k^ k"^ 

(3 - 2 + - ••• ) 



( 3 - 1 + 3 ^, - ... ) 

(3 - — + — 

'• 2 4! 



- ...) 



/. 00 ~ u(l + 0(k )) + v(l + 0(h )) 



(40) 

Q.E.D. 



This is to be expected, since the basis elements are the tensor product 
of two one- dimensional linear elements on rectangles. 
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A more interesting result is whether the same phenomenon happens if 
we use linear elements on triangles. We have the following result. The 
finite element approximation for (37) using the general linear element 
(j) = a + bx + cy on triangles is 



Ich f * * • • • • • 

12^^ Vi.j + ♦t.j+i ♦m.j + V.j-i + W.j+i 

'♦*+!, j+l ■ ♦jt.J+l' * ■ *11-1, J-l'* 

vie. 



This also has a fourth order phase speed accuracy as we now show. If we 

again look at 4>g.* as in (38) we get, for the first term in (41) 

^3 



21khoj 

12 



{3 + cos k + cos h + cos k cos h - sin k sin h} 



2ikhm r /• ,2 ,2 , , . hk^ . kh^ , i 

— — {6-k - h -kh+-jj- +'Ji' + ••• } 



while the second term is 



2ihu 



{3 sin k - sin h(l - cos k) - sin k(l - cos h)} 



3 

2ikhu r, ,2 ,2 ,,,hk, i 

—^-2 — {6-k - h -khH + . . . > 



and the third term is 



2ikv 



{3 sin h - sin k(l - cos h) - sin h(l - cos k)} 

. . 2ikhv r ^ .2 ,2 , , , k^h , 

{6-k - h - khH — h . . . J- 
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Thus 



o)*rw u(l + 0(hV)) + v(l + 0(hV)) (42) 

where m + n = 4. Thus we have shoxm 

Theorem 2: The finite element approximation to (37) using linear 

elements on triangles has fourth order accuracy in phase speed* 
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